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Nomenclature 

f melt fraction 

Ç specific heat at constant pressure 

H total volumetric enthalpy (J/m?) 

h average heat transfer coefficient (J/kg) 
k thermal conductivity (W/mK) 

L latent heat (J/kg) 

St Stefan number 

t time (s) 

T temperature (K) 

Tm melting temperature (K) 

Ts initial temperature of solid PCM (K) 
To constant temperature imposed on x=0 (K) 
P present nodal point 

ôt) position of liquid-solid interface 

X, y Cartesian coordinates 


1. Introduction 


Due to the attractive features of latent heat storage, phase 
change materials (PCMs) are mainly used to store energy at a fixed 
temperature (melting/solidification temperature) with high 
energy density |1]. PCMs have been applied in various systems 
and aspects such as energy storage systems, free ventilation, free 
heating and cooling for buildings, spacecraft, food, medicine 
conservation, smoothing the peaks of exothermic temperature in 
chemical reactions etc. PCMs are a better choice comparing to the 
sensible heat storage in applications, because of its nearly iso- 
thermal storing mechanism and high storage density. PCMs absorb 
abundant energy through the phase transition and release the 
stored energy. Table 1 compares the typical properties of different 
thermal storage materials tested at the room temperature. 
It indicates that PCMs can save 92.8% of mass and up to 90% space 
to store the same amount of thermal energy comparing to the 
sensible thermal materials such like concrete and water. 

Therefore, PCMs attract the researchers' interesting in practically 
individual and incorporation applications in different engineering 
fields. The possibility, feasibility, thermal performance and economic 


Table 1 
Properties comparison of different thermal storage materials. 


Property Unit Materials 
Concrete (sand Water Organic Inorganic 
and gravel) PCM PCM 
Density [kg/m3] 2240 1000 800 1600 
Specific heat [kJ/kg K] 11 4.2 2.0 2.0 
capacity 
Latent heat [kJ/kg K] - = 190 230 
Storage mass for [kg] 60,000 16,000 5300 4350 
10° kJ, avg 
Storage volume for [m3] 26.8 16 6.6 2:7 
10° kj, avg 
Relative storage 13.79 4 1.25 1.0 
mass 
Relative storage 10 6 2.5 1.0 
volume 


1. Storage mass and volume are calculated for storing 106 kJ energy with a 
temperature rise of 15 K for concrete (sand and gravel) and water. 
2. Relative mass and volume are based on the inorganic phase change materials. 


Greek symbols 

A difference 

oc thermal diffusivity of PCM (m?/s) 

n dimensionless number used in derivations as a tem- 
porary substitution 

A dimensionless number in solution to Neumann pro- 
blem, liquid fraction 

p density (kg/m?) 

Subscripts 

l liquid 

s solid 

i ith spatial step 

J jth time step 

eff effective 


analysis of using PCMs call a series of theoretical and experiential 
investigations, the mainly two methods used to research the thermo- 
physical parameters of PCMs and interrelated systems. 

The experimental approaches offer a better indication of the 
actual PCM behavior and performance in comparison to theore- 
tical analysis, as the experimental tests can present the PCMs' 
behaviors more directly, visibly and credibly without any pre-set 
assumptions however, the experiments are unachievable in some 
cases, such as the large scale or unsteady around environment, so 
there are still a few points need to be considered: 


è Time and cost consuming will be higher than a theoretical 
approach, hence, the budget and processes needs to be estab- 
lished and scheduled. 

© The scales of the experiment test rig have to be considered, and 
then suitable laboratory location and space needs to be 
determined to house the rig. 

© Test rig needs to be constructed and operated properly. 

© Proper environment parameters have to be simulated and 
adjusted to imitate the practical environment. 

© Relevant parameters need to be monitored, measuring appara- 
tuses need to be calibrated, and the failed data need to be 
eliminated. 


Except these agonizing experimental matters, there are still 
some unavoidable testing errors. However the theoretical methods 
can avoid all these weakness and predicate the PCMs performance 
conveniently. The major advantage of the theoretical/numerical 
approaches is that various conditions can be carried out by 
changing the variables in a numerical model. Therefore, more 
and more investigators prefer to study the phase change problems 
by mathematical solutions and numerical simulations. There are 
only eight review papers on PCMs that involves the aspect of 
mathematical solutions and/or numerical modeling of latent heat 
thermal energy storage (LHTES) have been published during 
2002-2012. Table 2 summaries these eight relative review papers 
and comments their mathematical and/or numerical aspects of 
LHTES. However, few logistic and comprehensive reviews on the 
mathematical solutions and numerical modeling for melting and 
solidification processes of PCMs were found in published litera- 
tures even though there are many individual publications in phase 
change problems, particularly in heat transfer mechanisms and 
simulations. Due to the complexity of phase transformations, 
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Table 2 
Some review papers related to numerical study of phase change problems. 


Reference Published Journal Comment 
year 
Zalba et al. [1] 2003 Applied thermal 
engineering 
Verma et al. [2] 2008 


energy reviews 


The review was divided into three parts: materials, heat transfer and applications. Heat transfer 

was considered mainly from a theoretical point of view, considering different simulation techniques. 
Renewable and sustainable Mathematical models of a LHTES were reviewed to optimize the material selection and to assist 

in the optimal designing of the systems. Some important characteristics of different models and their 


assumptions used were presented and discussed as well. 


Jegadheeswaran 2009 
and Pohekar [3] 
Zhu et al. [4] 2009 


energy reviews 
management 


Agyenim et al. [5] 2010 


Renewable and sustainable Various thermal conductivity enhancement techniques reviewed in this paper, and the issues related 
to mathematical modeling of LHTES with enhancement techniques are also discussed. 

Energy conversion and This paper presented an overview on dynamic characteristics and energy performance of buildings 
employing PCMs by three research methods used, i.e., simulation, experiment, combined 

simulation and experiment 

Renewable and sustainable This paper provided the formulation of the phase change problem. In terms of problem formulation, 


energy reviews it was concluded that the common approach has been the use of enthalpy formulation. 


Kuznik [6] 2011 

energy reviews 
Dutil et al. [7] 2011 

energy reviews 
Zhou [8] 2012 


Renewable and sustainable The review was the first comprehensive review of the integration of phase change materials in building 
walls. Various numerical studies concerning the integration were summarized. 

Renewable and sustainable The review presented models based on the first law and on the second law of thermodynamics. This paper 
tried to enable one to start his/her research with an exhaustive overview of the subject. 

Applied energy This paper reviewed previous works on latent thermal energy storage in building applications, including 


PCMs, current building applications and their thermal performance analyses, as well as numerical 
simulation of buildings with PCMs. 


a good understanding of the heat transfer mechanisms, phase 
change characteristics and the differences among various mathe- 
matical and numerical simulation methods is required before 
researchers starting their theoretical study. The aim of this paper 
is to comprehensively study the mathematical and simulation 
modeling applied for LHTES. Firstly it discusses the issues of the 
heat transfer mechanisms during the charging and discharging 
processes and the Stephan problem and Neumann problem. 
Secondly, it sums the strong and weak aspects of various mathe- 
matical solutions when conduction and convection heat transfer 
occurring inside the PCMs by the considering fixed grid method 
and the adaptive grid method. Finally, numerical modeling of 
various PCM packages’ geometries in terms of macro-encapsulated 
PCMs and micro-encapsulated PCMs are studied. 


2. Heat transfer mechanisms during the phase 
transformations 


2.1. Conduction and convection heat transfer 


During the charging and discharging processes, the possible 
heat transfer mechanisms are conduction and convection. How- 
ever, the issue of which heat transfers mechanism takes the main 
role in different stages of phase transformation has been argued 
for decades. When PCMs are used to store or release thermal 
energy, conduction is usually believed playing the most important 
role on the heat transfer during the phase transformation pro- 
cesses [2-4]. As early as 1831, Lamé and Clapeyron have conducted 
the first study on phase change problem by only considering the 
pure conduction. However, some researchers persist that natural 
convection is a more important mechanism in the phase change 
process especially in the melt region. Lazaridis [100] proposed a 
study of the relative importance of conduction and convection in 
1970. A pioneer study performed by Sparrow et al. in 1977 [41], 
they concluded in their study that natural convection could not 
be ignored in the analysis of phase change problems. In 1994, 
Hasan [5] concluded that the convection heat transfer active an 
important role in the melting process, and a simplified model by 
only considering the conduction heat transfer does not describe 
the melting process properly. Later, Lacroix et al. [6] reported 
the similar findings in their research that natural convection is 
the main heat transfer mechanism during the melting process. 


In 1999, Velraj [7] obtained the same conclusion in their research 
and reported that during the melting process natural convection 
occurs in the melt layer which results in the heat transfer rate 
increase comparing to the solidification process. Buddhi et al. [8] 
proposed an explanation for this phenomenon that the density 
differences between the solid and liquid PCM induce the buoy- 
ancy, which causes the convective motions in the liquid phase. 
However, Zhang and Yi [9] believed that with the solid PCM 
melting into liquid, the PCM volume keeps increasing, which 
results in the convection becoming the predominant heat transfer 
mechanism rather than conduction. Sari et al. [10] found that the 
heat transferred from a heat exchanger to a PCM of stearic acid is 
largely influenced by natural convection at the melting layer, in 
addition to conduction and forced convection heat transfer. 

For the melting process, the PCM changes its phase from solid 
to a mushy state, and then liquid, which is reversible during the 
solidification process. Hence there are six stages to finish a 
charging and discharging cycle. Therefore it is possible in certain 
stages of the phase transformation process there are more than 
one kind of heat transfer mechanisms acting at work, but how to 
weight the percentages of conduction and convection heat transfer 
in each stage has been the main challenge for the researchers 
currently. This paper will review the most common research 
methods used for the conduction and convection heat transfer 
inside the PCM packages. 


2.2. Stefan problem 


Moving boundary problem named Stefan problem is another 
issue to develop a numerical modeling of PCMs [11,12]. The 
simplest of the Stefan problems is the one-phase Stefan problem 
since only one-phase involved. The term of ‘one-phase’ designates 
only the liquid phases active in the transformation and the solid 
phase stay at its melting temperature. Stefan's solution with 
constant thermophysical properties shows that the rate of melting 
or solidification in a semi-infinite region is governed by a dimen- 
sionless number, known as the Stefan number (St), 


[CTi = Tn) 


St L 


(1) 
where C; is the heat capacity of the liquid PCMs, L is the latent 
heat of fusion, and T; and Tm are the surrounding and melting 
temperatures, respectively. 
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The most published solutions are to solve the one-dimensional 
cases subjected to an infinite or semi-infinite region with simple 
initial and boundary conditions. Whereas with the heat transfer 
continuing, the interface boundary is constantly moving as the 
liquid and solid phases shrinking and growing, which disable the 
prediction of the boundary location | 13]. Because of that the solid- 
liquid interface is not fixed, but moving with time, the heat 
transfer mechanisms during a PCM phase transformation process 
are complex. Therefore the phase change transition is difficult to 
analyze owing to the three reasons: the solid-liquid interface is 
moving; the interface location is non linear; it consists of thermal 
conduction and natural convection heat transfer mechanisms. 
Since these three factors, the non-linearity of the governing 
equations is introduced to the moving boundary, and the precise 
analytical solutions are only possible for a limited number of 
scenarios [14]. This view is shared by Kürklü et al. [15] who 
pointed out that the mathematical models have been proposed in 
each research only applied to very specific boundary conditions, 
hence are not feasible to more complex practical applications. This 
Stefan problem is further complicated by the fact that most of the 
methods previously proposed by researchers only involve one 
moving boundary, whereas it actually consist of more than one 
interface location [16]. 


2.3. Neumann problem 


The Stefan problem was extended to the two-phase problem, 
the so-called Neumann problem which is more realistic [17]. 
In Neumann problem, the initial state of the PCM is assumed to 
be solid, during the melting process, its initial temperature does 
not equal to the phase change temperature, and the melting 
temperature does not maintain at a constant value. If the melting 
happens in a semi-infinite slab (0O<x<oo), the solid PCM is 
initially at a uniform temperature (T;<Tm), and a constant 
temperature (Tg) is imposed on the slab surface x= 0, with the 
assumptions of constant thermo-physical properties of the PCM, 
the problem can be mathematically expressed as follows: 

Heat conduction in solid or liquid region 


pC—= ka, (2) 


where p is the density, C is the specific heat, k is the thermal 
conductivity, and t and x are the time and space coordinates 
respectively. 

The heat fluxes transferring from the liquid phase to the solid- 
liquid interface, as well as the latent heat absorbing rate by the 
melting PCM, the movement of the solid/liquid interface can be 
determined through the following energy balance: 


OXt aT; Ts 
pi = kı tk 


JK (3) 


where L is the latent heat of fusion of the PCM and p is the density 
of liquid PCM. 
Initial condition 


T(x > 0,t =0) =T; <Tm (4) 
Solid-liquid interface temperature 

T(x = 6(t), t>0)=Tm (5) 
With the following boundary conditions: 

TO, t)=To>Tm for t>0 (6) 

T(x, =T; for x>œ, t>0 (7) 


where 6(t) is the position of the solid—liquid interface (melting 
front). Fig. 1 clearly illustrates this two-phase Stefan problem. 


T(0,t) = To > Tm 


Fig. 1. Schematic illustration of the two-phase Stefan problem. 


Analytical solution to such a problem was obtained by Neu- 
mann in term of a similarity variable 7 


ôt) 
2/ C] 
The final Neumann's solution can be written as 
Interface position 


&(t) = 2A,/ «it (9) 


Temperature of the liquid phase 


q= (8) 


erf(x/2/ « ıt) 
T(x, t) =T,- (Tı—Tm) erf (10) 
Temperature of the solid phase 
insti eee an 


erfa a1] &s) 


The A in Eq. (9)-(11) is the solution to the transcendental 
equation 
A y Cs = AVT 
exp erf(A) SEI exp oÀ" / o sjerfeAy x 1/ xs) 


a2) 
where 
St; = CiTi—-Tm) (13) 
L 
St; = — Ts) (14) 


However, the Neumann's solution is applicable only for moving 
boundary problems in the rectangular coordinate system. 


2.4. Other possibility mechanisms in phase change process 


Furthermore, there are several other issues with the use of a 
theoretical approach in the study of PCMs. Alexiades [18] pointed 
out that there were many mechanisms involved during a PCM 
phase transition, such like a change in volume, density, thermal 
conductivity, specific heat capacity, super-cooling, etc. Conse- 
quently, accurate reflection of the proposed mathematical solution 
and numerical model need to consider the dynamic properties of 
the phase change process. Another major issue with PCMs is that 
they act as self-insulating materials. When PCM solidification 
occurs from the top of the heat surface, solid insulating layer will 
be developed which moves inward during the whole solidification 
process. With the increase in the size and thickness of the solid 
layer, the heat transfer rate from the liquid PCM to the heat 
exchanger surface decreases until it becomes so small that will not 
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be possible to maintain at an acceptable heat transfer rate. The 
self-insulating characteristic of the PCM becomes a strong problem 
when bulk containerization is used. For example the large cylind- 
rical used as containers filled with PCM, the self-insulating issue 
can become very serious. Large heat exchange area will accelerate 
the formation of the crust and heavily decrease the heat transfer 
rate. Radhakrishnan et al. [14] analyzed that the method to lighten 
the crust issue is minimizing the depth of the PCM pack. So that 
even a large majority of the PCM performs solidification, the 
thermal resistance is reduced to allow the removal of the latent 
heat released from the solidified PCM. The self-insulating char- 
acteristic of the PCM during discharging process needs to be 
considered in any mathematical model as this can affect thermal 
energy discharging time and rate. Furthermore, the self-insulating 
has a direct effect on the phase transition boundary, where the 
solid and liquid region meets and co-exists, which is linked back to 
the Stephan Problem and Neumann's method again. 

Logical heat transfer mechanisms during the charging and 
discharging periods of the PCM are the essential elements to 
develop any accurate mathematical model. However, as described 
there are quite a few issues to properly describe the heat transfer 
mechanisms: conduction and natural convention during the phase 
change process. This paper reviews most of the researches 
employing the mathematical calculation and numerical simulation 
published in the articles and reports, but does not find any 
consensuses on which one is the dominant heat transfer mechan- 
ism, which one can be ignored and how to combine these two 
transfer mechanisms in each modeling. Different researchers have 
developed their specific methods to evaluate the thermal energy 
transfer during the phase transition process for various situations 
[19-23]. 


3. Mathematical solutions 


The analyzes of the heat transfer mechanisms in phase change 
processes are complex owing to the movable solid-liquid bound- 
ary, which depends on the latent heat absorbing or releasing rate 
at the interface. In order to simplify the simulation of the phase 
change problems, assumptions have been made by various 
researchers in their studies. In this section, different numerical 
simulations in terms of conduction heat transfer and conduction 
and/or convection heat transfer are discussed and analyzed. 


3.1. Conduction acting as the only heat transfer mechanism 


To develop a mathematical model, the dominant heat transfer 
mechanism whether conduction or convection, during the char- 
ging and discharging processes, need to be analyzed and deter- 
mined. Although there are quite a few issues in the heat transfer 
mechanisms during the phase transformation, some researchers 
simplify this process by only considering conduction in pure sub- 
stances. In this case, the mathematical solutions can be classified 
as follows: fixed grid and adaptive grid. 


3.1.1. Fixed grid 

Fixed grid method does not require explicit treatment of 
conditions on the phase change boundary, and just requires the 
use of fixed-grid schemes able to cope with strong nonlinearities. 
Therefore this method is able to utilize standard solution proce- 
dures for the energy equations. They are amenable to physical 
interpretation and easy to implement, especially for multidimen- 
sional problems and for mufti-interface problems. 

In the fixed grid method, the heat flow equation is approximated 
by finite difference replacements for the derivatives in order to 
calculate the values of temperature T;j, at any time tn =nAt, the 


Fig. 2. Position of the moving boundary in a fixed grid. 


moving boundary will be located between two adjacent grid points, 
for instance, between mAx and (m+1)Ax (at a location of eAx, 
where 0 <e< 1), as illustrated in Fig. 2. 

The mathematical solution of the one-phase ice-melting pre- 
sents a simple illustration of this method. The instantaneous 
temperature at point (i,j+1) is calculated from the temperatures 
of the previous step on the basis of the following formulation: 


At : 
Ty = Tut (Ea) T- 2Tij+Ti+1j) i=0,1,...m—1. (15) 


In terms of three-point Lagrange interpolation instead of 
Eq. (14), the temperature at x = mAx is computed as, 


2At 1 1 
T =Tnnt+ Lms T ) 16 
mn+1 m,n (x) ( m—1,n 5 m,n ( ) 


En+1 Er 


The variation of the location of the moving boundary is 


At E 1 
En+1 = En (xe) (Enin Tnn) (17) 


As shown, the numerical solution of this fixed grid method is 
carried out on a space grid remaining at a fixed value throughout 
the mathematical calculation. To approximate the Stefan condi- 
tions both for the moving boundary and the partial differential 
equation at the adjacent grid points, various mathematical solu- 
tions have been proposed. Murray and Landis [24] proposed two 
fictitious temperatures: one is achieved by quadratic extrapolation 
from the temperatures in the solid PCM and the other one is from 
the temperatures in the liquid PCM. The melting temperature and 
the position of the moving interface are incorporated in the 
fictitious temperatures, which are then substituted into a stan- 
dard approximate to calculate the temperature near the inter- 
face instead of using special formulae. Ciment and Guenther [25] 
introduced a method of spatial mesh refinement on both sides of 
the moving boundary. With this method, Lazaridis [26] applied the 
explicit finite difference approximations to solve two-phase soli- 
dification problems in both two and three space dimensions. 

The major advantage of the fixed grid method is that the 
mathematical calculation of the phase transition can be achieved 
through simple modifications of existing heat transfer numerical 
methods and/or software. As such, they can be used for modeling of 
a variety of complex moving boundary problems. Basu and Date [27] 
and Voller et al. [28] reviewed the applications of the fixed grid 
methods for phase change problem. The two mainly fixed-grid 
methods that have been proposed are enthalpy-based methods 
and the effective heat capacity method. However, the fixed grid 
method sometimes cannot work efficiently when the boundary 
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h 


phase change 
(melting) 


0 Tn T 


Fig. 3. Relationships between enthalpy and temperature for isothermal phase 
change. 


moves along a distance larger than a space increment in a time step, 
depending on the velocity of the moving boundary. Therefore, the 
adaptive grid is developed by researchers and will be discussed in 
Section 3.1.2. 


3.1.1.1. Enthalpy based method. The existence of nonlinearity of the 
liquid-solid interface and the unknown location of the moving 
boundary are the two main challenges in simulating the phase 
change problems [29]. To deal with these, a new model formu- 
lation called the enthalpy method was introduced [30]. In the 
enthalpy models, the enthalpy is used as a dependent variable 
value along with temperature. By introducing an enthalpy method, 
the phase-change problem becomes much easier since the 
governing equations are the same for the two phases. The 
interface conditions are automatically achieved and a mushy 
zone between the two phases is created. This zone avoids sharp 
discontinuities that may create some numerical instability. Hence, 
the enthalpy mathematical methods are most attractive and 
commonly used to handle phase change problems where a 
solid-liquid mush zone is present between two phases. Hunter 
[31] confirmed that the enthalpy method is most suitable for 
typical applications, the reason for that is this method does not 
require explicit treatment of the conditions on the phase change 
boundary [32]. 

Enthalpy-based method can be illustrated by considering a 
one-dimensional heat conduction controlled phase transforma- 
tion. For a phase change process, energy conservation can be 
expressed in terms of total volumetric enthalpy and temperature. 
This relationship between the enthalpy and temperature is usually 
assumed to be a step function for isothermal phase change problems. 
Fig. 3 shows the enthalpy-temperature curves for isothermal phase 
change. 

For isothermal phase change, the conservation of energy can be 
expressed in terms of temperature and the total volumetric 
enthalpy [33] 
oH 
a V[k(VT)] (18) 
where the total volumetric enthalpy H is the sum of sensible and 
latent heat, and can be expressed in terms of temperature of the 
PCM as follows: 


HT) = AHT) +p, fA a9 


The first term on the right side of Eq. (19) accounts for the sensible 
heat of the PCM while the second term accounts for the latent heat 
of the PCM. 


And where 
T 
AH = pcdT (20) 
Tm 


To be able to calculate the total enthalpy, the liquid fraction J is 
given as follows: 
0 T<Tm_ solid 
A= ]0,1[ T=Tm mushy (21) 
1 T>Tm liquid 


Integrating the Eqs. (18) and (20), an alternate form for one- 
dimensional heat transfer in the PCM was obtained: 


aAH 2 (aP AHN i jf (22) 


ot Ox Ox ot 


where æ is the thermal diffusivity. 


3.1.1.2. Effective heat capacity method (energy based method). The 
effective heat capacity method is also a common method used in 
modeling phase change problems except for the enthalpy method. 
The effective heat capacity method was first proposed in [34,35]. In 
the effective heat capacity method, the latent heat is approximated 
by a large heat insensible form over the phase change temperature 
interval, (T2 — T1) [36]. The effective heat capacity of the PCM (Cerf) 
is directly proportional to the energy stored and released during the 
phase change but inversely proportional to the interval of the 
melting or solidification temperature range. During the phase 
change the heat capacity of the PCM is 


L 
Ceff =p- "F (23) 


where T, is the temperature at which melting or solidification 
begins and Tz is the temperature at which the PCM is totally melted 
or solidified. The following is a definition of the effective heat 
capacity for each phase change period. 


Cs T<T; solid 
Cop =} r rts Ta <T<T2 mushy 24) 
Ci T>T2 liquid 


In terms of the definition of the effective heat capacity, the 
energy equation for one dimension can be written as follows: 


or æT 
Plott ae = koa (25) 


Although the enthalpy based method and the effective heat 
capacity method achieve much simpler mathematical operations 
with reasonable accuracy, it is worth noting that each of the 
methods has its inherent disadvantages. Enthalpy method is 
difficult to handle super-cooling and temperature oscillation 
problems. Whilst the effective heat capacity method is difficult 
in resolving the phase change problems when the phase change 
temperature range is small, and is not applicable for the cases 
where phase change occurs at fixed temperature. 

Two approaches, finite-difference and finite-element techniques, 
are firstly used to solve the phase change problems numerically. 
The first finite-difference studies were carried out by Lamé and 
Clapeyron in 1831 and Stefan in 1891, concerning ice formation. 
In 1912 the results of Neuman were published, establishing the precise 
solutions to more general phase change problems. In 1943, London 
and Seban [37] analyzed the process of ice formation for different 
geometries (cylinder, sphere, and flat plate). However, Shamsundar 
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and Srinivasan [38] asserted that London's one-dimensional formula- 
tion led to errors by increasing the progress of the solidification 
process and proposed a two-dimensional formulation for cylinders. A 
one-dimensional formulation was presented by Bonacina et al. to 
assess the accuracy of the effective heat capacity method [39]. Rabin 
and Korin presented a multidimensional numerical solution based on 
the effective heat capacity method for solving transient phase change 
problems by using the finite-difference method [40]. The proposed 
mathematical method showed a good agreement with two exact 
solutions and two numerical solutions based on finite-difference and 
finite-element methods. 

Although the finite-difference method has been traditionally 
employed for solving phase change problems in the past, oscilla- 
tions in the front position and/or temperature occur due to a poor 
heat balance when the phase change front lies between nodes ina 
finite difference solution [41,42]. Nowadays the finite-element 
method has been more attractive, because of the ability of the 
finite-element method to handle complex coupled thermomechanical 
mechanisms with various and complex boundary conditions [43]. 
Consequently, finite-element methods can either accurately repre- 
sent the temperature history or track the phase front, or both. A 
few earlier studies based on finite-element formulations were 
carried out by Rubinsky and Cravahlo in 1981 [44] and Ettouney 
and Brown in 1982 [45], concerning on the location of the solid- 
liquid interface. To deal with oscillations, Comini et al. [46] and 
McAdie et al. [47] developed finite-element based enthalpy 
formulations, and applied the three level Predictor—Corrector 
and Regula-Falsi methods to account for the sudden jumps or 
sharp changes in the enthalpy-temperature relationships. Bhatta- 
charya et al. developed a new enthalpy based method to study 
phase change in a multicomponent material using the Galerkin 
finite-element method. This method can efficiently capture multi- 
ple thawing fronts which may originate at any spatial location 
within the sample [48]. Hence, finite-element methods are pre- 
ferred to be used in enthalpy approaches and the effective heat 
capacity method. 

Beside the finite-difference and finite-element methods, the 
finite-volume and control-volume- finite-element methods were 
also employed to handle the phase change problems [39]. 


3.1.2. Adaptive grid 

The problem associated with the fixed grid method can be deal 
with by utilizing the adaptive grid method to improve the 
computational efficiency, using the so-called adaptive grid (or 
front tracking) numerical schemes to capture the moving bound- 
ary [49]. In the adaptive grid method, the exact location of the 
moving boundary is evaluated on a grid at each step. There are two 
main approaches used to achieve this: the interface-fitting grids 
and the variable-space grids. The interface-fitting grids (also 
referred to as the variable time step methods), where a uniform 
spatial grid but a non-uniform time step are used, has been 
repeatedly employed to solve two-phase and one-dimensional 
problems. The another approach is variable-space grids, also 
known as the dynamic grids, where the number of spatial intervals 
is kept constant and the spatial intervals are adjusted in such a 
manner so that the moving boundary lies on a particular grid 
point. Thus, in these methods the spatial intervals are a function of 
time. The applications of the variable grid method to study the 
phase change problem can be found in [50-54]. 

Lacroix and Voller [55] performed a study on simulation of 
diffusion/convection controlled solidification processes in a rec- 
tangular cavity by using the finite-difference techniques. They 
concluded that the fixed grid must be finer for solving the melting 
problem of a material with a unique melting temperature. Com- 
parison between the fixed and adaptive grid has also been done by 
Bertrand et al. [56]. They found that the adaptive grid method is 


better adapted to the melting problem than the fixed grid method. 
However the adaptive grid method may fail to simulate the 
situations where the transition from the liquid to the solid phase 
is not a macroscopic surface. 

Table 3 presents some references of the mathematical and 
numerical studies undertaken with the only consideration of the 
conduction heat transfer. 


3.2. Conduction and convection heat transfer mechanisms 


During the phase transformation especially melting process, 
the temperature and concentration gradients in the liquid phase of 
PCM keep varying, which results in the movement of the liquid 
PCM, named convection occurring under the action of buoyancy 
forces due to the density gradients. In the previous numerical 
investigations, the convection heat flow in the liquid phase 
received less attention than conduction owing to the limited 
computational capabilities of computer and the mathematical 
complexities to formulate the convection heat transfer during 
the phase transformation. The influence of the natural convection 
on the phase change problems were initially considered in 1977. 
The study of the possible role of the natural convection in the melt 
region has been conducted by Sparrow et al. in 1978 [67]. The 
findings of this study indicated that the natural convection is first- 
order important and has to be accounted in the phase transforma- 
tion analysis. A number of researchers [68-71] have also reported 
that the convection affects not only the rate of melting or 
solidification but also the resulting of the distribution of the liquid 
phase of the PCM in a multi-component system. One of the 
pioneer numerical studies was carried out by Yao and Chen in 
1980 [72], by using an approximate solution. They studied the 
effect of the natural convection on the phase change process and 
concluded that it strongly depends on the Rayleigh number. The 
equation of ke/k,;=cRa" includes an effective of thermal conduc- 
tivity was used to describe the influence of natural convection on 
the melting process [73,74]. 

The influence of the convection heat transfer on the melting 
and freezing processes has been studied by Lamberg et al. [64], 
they found that when the effect of natural convection is omitted in 
the simulation. The numerical result was poor compared with the 
experimental result during the melting process, whilst it showed a 
good estimation during the freezing process. Hence it is essential 
to model natural convection in the liquid PCM during the melting 
process. Besides the numerical methods mentioned in Section 3.1, 
the integral method [75], boundary fixing method [76], unstruc- 
tured finite-element method [77], enthalpy-porosity approach 
[78-80], coordinate transformation method [81,82] and equivalent 
thermal capacity method [83] have been developed for natural 
convection simulation. 

Yao and Cherney [84] studied the effect of the natural convec- 
tion on the melting of a solid PCM around a hot horizontal cylinder 
by using the integral method. The results demonstrated that the 
integral solution had surprising accuracy when it was compared 
with the quasi-steady solution when Stefan number was small. 
Rieger et al. [85] investigated numerically the melting process of a 
PCM inside a heated horizontal circular cylinder. Both heat con- 
duction and convection have been taken into account to treat this 
moving boundary problem, and the complex structure of the time- 
wise changing physical domain (melt region) have been success- 
fully overcome by applying body fitted coordinate technique. 
Ismail and Silva [86] developed a 2D mathematical model to study 
the melting of a PCM around a horizontal circular cylinder 
considering the presence of the natural convection in the melt 
region. A coordinate transformation technique was used to fix the 
moving front. The numerical predictions were compared with 
available results to establish the validity of the model, and a 
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Table 3 
Numerical studies by only considering the conduction heat transfer. 


Reference(s) Dimensional Numerical Assumptions Comment Validation 
simulation method 


Silva et al. 1D Enthalpy based e 1-D conduction controlled the phase The model can be used to predict the dynamic Yes 
[57] method change process. performance of this kind of vertical rectangular 
e The fluid flow was fully developed with heat exchange unit with reasonable accuracy 
temperature independent properties. 
e Thermal radiation between the wall 
and fluid was neglected. 
Kürklü et al. 2D Energy based e The phases were homogeneous. Calculating the temperature of air and PCM Yes 
[15] equation e Conduction heat transfer was at the formerly defined volume nodes. 
only considered. 
e Heat loss or gain from the store was 
neglected. 
Costa et al. 2D Energy based e The thermophysical properties of the PCM Calculations have been made for the melt fraction Numerical 


[3] equation. and fin were independent of temperature. and energy is stored by conduction only [33,58] 
Fully implicit The PCM was initially in the solid phase. 


finite-difference e The PCM was homogenous and isotropic. 
method. 
Gong et al. 2D Enthalpy based e The heat-transfer fluid was incompressible For mode1, simulation has been done by No 
[59] method. and viscous dissipation was negligible. introducing 
Standard Galerkin e The fluid flow was radially uniform and hot and cold fluids from the same end of the tube. 
finite-element the axial velocity was an independent While for mode 2, simulation has been done by 
parameter. introducing hot and cold fluids from the different 
e Thermal losses through the outer wall ends of the tube. 


of the PCM were negligible. 
Conduction heat transfer was the 
mainly transfer method. 


Sharma et al. 2D Enthalpy based e The thermo-physical properties of PCM's and Calculations were made to study the effect of No 
[60] method fin material were independent of thermal conductivity and thermal capacity of fatty 
Fully implicit temperature but different for solid and liquid acids, and conductivity of heat exchanger materials 
finite-difference phases. and their effect on melt fraction. 
method e The PCM was initially in solid phase. The 
PCM was homogeneous and isotropic. 
e The mode of heat transfer was 
conduction only. 
Dwarka and 3D Implicit enthalpy e Air temperature was constant. Comparison of the thermal performance of No 
Kim [61] based method. e Initial temperatures were the same through randomly mixed and laminated PCM drywall 
the board. system 
e There was no energy loss to surroundings. 
Conduction heat transfer was only 
considered. 
e All thermo-physical properties were 
constant except for the heat capacity. 
Dolado et al. 1D Finite-difference e Only conduction heat transfer was A model was developed to simulate the Yes 
[62] method considered inside the PCM plate. performance 
e The temperature of the airflow was analyzed of a LHTES, and analyze the heat transfer between 
in a one dimensional way. the air and a commercial macroencapsulated PCM. 
Fukai et al. 2D Enthalpy based e The surfaces of the heat exchanger were A numerical model predicted well the experimental Yes 
[63] method adiabatic. outlet fluid temperatures and the local 
Control-volume e The cross-sections of the tubes was square temperatures 
method e The convective heat transfer term was in the composite. 
neglected. 
Lamberg et al. 2D Enthalpy based e The heat conductivity and density of the PCM Both numerical methods gave good estimations for Yes 
[64] method and Effective and container were constant. the temperature distribution. Particularly, when 
heat capacity e The physical properties chosen for the PCM the temperature range was 2 °C, the effective heat 
method were the average values of solid and capacity method was the most precise numerical 
liquid PCM. method. 
e The convection heat transfer was negligible 


during solidification process. 


Zivkovic and 1D Enthalpy based e The effects of natural convection within the The cylindrical container required nearly twice of Yes 
Fujii [65] method melt were negligible and can be ignored. the melting time as for the rectangular container of 
e The PCM behaved ideally. the same volume and heat transfer area 
e The PCM was assumed to have a definite 


melting temperature. 
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Table 3 (continued ) 


Reference(s) Dimensional Numerical Assumptions Comment Validation 
simulation method 
Saman etal. 2D Enthalpy based e The thermo physical properties of PCM The numerical model was able to predict quite Yes 
[66] equation can be assumed constant. accurately the melting time and the heat transfer 


The PCM was initially solid and its 


rate during the melting with the geometry being 


temperature was assumed at a certain value considered 


below the melting point. 


The natural convection was taken 


into account during charging period. 


satisfactory agreement was found. They concluded that the 
numerical model was adequate to represent the physical situation 
of the proposed system. 

Furzerland [87] investigated the enthalpy method and the 
coordinate transformation method through the comparison of 
the solutions of specific problems of one dimensional heat transfer 
by considering pure convection. One of his conclusions was that 
the enthalpy method is very attractive owing to these: it is easy to 
program and there are no computational overheads associated 
with tracking the moving interface within a specific range of 
fusion temperatures. 

A summary for various numerical solutions on phase change 
problem considering both conduction and convection is presented 
in Table 4. 


4. Numerical models for various PCMs capsules and packages 


It can be easily found that solid-liquid PCMs have been widely 
used for studies on phase change problems compared to solid- 
solid, liquid-gas and solid-gas PCMs. This is due to their large heat 
storage capacity and little volume change during the phase change 
process. Since solid-liquid PCMs experience a phase transition 
during the charging process and/or charging process, encapsula- 
tion is required for holding the liquid and/or solid phases of the 
PCM and keeping it isolate from the surrounding. This ensures the 
flexibility in frequent phase change processes, and an increase in 
heat transfer rate. It is worth noting that the heat transfer 
characteristics in the encapsulated PCM would change signifi- 
cantly depending on the parameters of various encapsulations. 
Hence, vast of numerical studies on phase change problems have 
extended to engineering applications. In this section, a compre- 
hensive review on numerical models for encapsulated PCMs in 
terms of macro-encapsulated PCM and microencapsulated PCM is 
presented. 


4.1. Macro-encapsulated PCMs numerical models 


Macro-encapsulation is a common form of encapsulation used 
for thermal storage applications. The shape of macro-encapsulate 
varies from rectangular panels, spheres to cylinders. 


4.1.1. Rectangular encapsulation 

The most intensely investigated LHS containment is the rec- 
tangular encapsulation because to its simple boundary conditions 
and easy work principle. The first numerical study on rectangular 
geometry has been conducted by Shamsundar and Sparrow 
[107,108] in 1975 and 1976 by applying the finite-difference 
method to investigate the solidification of a flat plate. Saamsundar 
also used the same method to the case of a square geometry [109] 
in 1978. 

Hamdan and Elwerr [110] presented a simple 2-D numerical 
model to study the melting process of a solid PCM within a 


rectangular enclosure. In this model, the rate of melting depends 
essentially on the properties of the PCM, such as the thermal 
diffusivity, viscosity, conductivity, latent heat of fusion and specific 
heat. Natural convection was treated as the dominant heat transfer 
mode within the melt region. Conduction was only assumed to 
take place within the layer very close to the solid boundary. The 
numerical results showed a high agreement with previously 
published experimental results [111,112]. Therefore, it can be 
effectively used to predict the energy storage performance of the 
PCMs contained within the rectangular encapsulation. 

Further research was carried out by Lacroix [113] in 2001, who 
presented a mathematical model based on the energy conserva- 
tion equation to study the contact melting of a sub-cooled PCM 
inside a rectangular cavity including the natural convection effect. 
Numerical results demonstrated that the melted fraction at the 
bottom of the cavity is larger, by an order of magnitude than that 
of the conduction dominated melting at the top. The melting 
process is essentially governed by the magnitude of the Stefan 
number and strongly influenced by the lateral dimensions of the 
cavity. Jiji and Gaye [114] also analytically examined one- 
dimensional solidification and melting of a slab with uniform 
volumetric energy generation. They found that the low thermal 
conductivity of the PCMs present a significant challenge in the 
design of PCM-based electronic cooling systems. 

In order to address the inherent drawback of the low thermal 
conductivity of PCM, various numerical solutions have been 
applied for PCM with different thermal conductivity enhancers 
(TCEs). Lacroix and Benmadda [115] numerically analyzed the 
melting of PCMs in a rectangular enclosure with horizontal fins. 
A 2-D enthalpy model was used to solve the phase change 
problem using a fixed-gird method. The effects of the number 
and length of fins on the melting rate were considered in this 
study. It was concluded that a few longer fins were more effective 
in the melting rate than a large number of shorter fins. The effect 
of vertical fins on melting rate was studied by the same team [116] 
by using a fixed-grid enthalpy model in 1998. It was found that the 
onset of natural convection in the melted zone was delayed when 
the distance between the fins was decreased. According to their 
results, the optimum fin spacing decreases as the Rayleigh number 
increases. Akhilesh et al. [117] numerically studied a rectangular 
module with vertical fins, which was heated from above. The 
numerical model was solved by using finite-difference method. 
Only conduction heat transfer was included in this study. The 
analysis showed that more fins increase the rate of energy stored 
in the melting PCM until a critical value. There is no obvious 
performance enhancement by only increasing the number of fins 
beyond the value. Similar results were obtained in another 
research of Gharebagi and Sezai [118] in 2005 considering an 
enclosure with vertical fins added to a horizontal heated wall. The 
results indicated that the heat transfer rate to the melting PCM can 
be improved by adding fins. In addition, vertically heated walls 
with horizontal fins exhibit better performance than horizontally 
heated walls with vertical fins. 
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Numerical studies by considering both the conduction and convection. 


Reference(s) Dimensional Numerical solution method 


Costa et al. 1D 
[3] 


Halawa et al. 2D 
[88] 


Wang etal. 2D 
[90] 


Jana [91] 2D 


Hamdan and 3D 
Al-Hinti 
[93] 


Zhao et al. 1D 
[95] 


Liu et al. [97] 2D 
3D 


Shmuelietal. 2D 
[99] 


Energy based equation 
Fully implicit finite- 
difference method 


Enthalpy based method 


Energy based method 
Finite-volume method 


Enthalpy-porosity based 
method 
Finite-volume method 


Energy based equation 


Enthalpy porosity method 


Enthalpy porosity method 


Assumptions 


The thermophysical properties of the 
PCM and fin material were independent 


of temperature. 


The PCM was initially in the solid phase. 
The PCM was homogenous and isotropic. 


For the PCM melting process, the PCM 
was solid and its temperature was 
assumed at a certain value below the 


melting point. 


For freezing process, the PCM was 
initially liquid and its temperature was 
assumed at certain value above the 


melting point. 


The PCM was pure and homogeneous. 
The melting front was an isothermal 


interface in reality. 


The thermophysical properties of the 
PCM were constant, with the exception 
of the linear density-temperature 


relation. 


The natural convection flow of the liquid 
phase was incompressible and laminar 
with no viscous dissipation. 

The density of PCM was constant in 


this study. 


The density was considered to be 
constant in the unsteady and convective 
terms and was allowed to vary only in 


the body-force term, 


The phase change occurred at a single 


temperature, 


The temperatures of liquid and solid at 
the interface were equal. 


The solid-liquid interface was 


smooth, plane. 


The solid phase was initially at the 


melting temperature. 


The effect of inertia inside the thermal 
boundary layers was negligible relative 
to that of buoyancy and viscosity. 

The heat loss from the walls of the 
enclosure was negligible. 

The PCM material was pure, homogenous 


and isotropic. 


Boussinesq approximation was used in 


this analysis. 


The heat transfer process inside the 
molten layer was one-dimensional, and 
the tangential temperature gradient was 


neglected. 


The thermo-physical properties of PCM 


were constant. 


The process was quasi-steady, and the 
acting forces on PCM were balanced at 
every moment of the melting. 


The two-phase mixed region can be 
described with the porosity. 

The effective thermal conductivity of the 
mixture can be calculated as the volume 
average of the conductivities of porous 
matrix material and PCM 

The porous media and fluid were not 
assumed to be in thermal equilibrium. 


Both solid and liquid phases were 
homogeneous and isotropic. 
The melting process was axisymmetric. 


Comment Validation 
Calculations have been made for the melt Numerical 
fraction and energy stored under conduction [30,58] 
plus convection condition. 


Typical characteristics of the melting/freezing Numerical [89] 
of PCM slabs were discussed. Considerations 
in the design of the LHTES were also given. 


Calculations were performed to study the No 
effect of thermal conductivity and thermal 
capacity of fatty acids and conductivity of 

heat exchanger materials and their effect 

on melt fraction. 


Prediction of solidification and melting Numerical and 
processes of pure materials Experimental 
using moving grids. [92] 


The melting process of a solid PCM contained Experimental 
in a rectangular enclosure heated from [94] 

a vertical side at a constant heat flux 

was investigated analytically. 


Numerical and 
Experimental 
[96] 


The expression of melting velocity was 
obtained by using the lubrication theory. 


The effects of the structural parameters of the Experimental 
porous media and inlet conditions of HTF on [98] 

the thermal performances of LHTES were 

numerically analyzed. 


The findings of the present study made it Experimental 
possible to define the heat transfer [100] 
mechanisms from the tube wall to the liquid 
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Reference(s) Dimensional Numerical solution method Assumptions 


The molten PCM and the air were 
incompressible Newtonian fluids, 


Comment Validation 


PCM and then to the solid phase at various 
locations and instants. 


and laminar flow was assumed in both. 


Ye et al. 2D Enthalpy porosity e Density and dynamic viscosity of the The numerical study provided a detailed Experimental 
[101] approach PCM depended on temperature. knowledge regarding interface heat transfer [102] and 
e Variation of property values in liquid and rate provides a deeper understanding the Numerical 
solid phase PCM was imposed using heat [103] 


piecewise linear interpolation. 
Constant thermophysical properties 


transfer mechanisms. 


were specified for aluminum. 


Tao and He 2D 
[104] 


Enthalpy method 


axisymmetric model. 


The axial heat conduction and viscous 
dissipation in the HTF was negligible. 
The flow of HTF was treated as one 
dimensional fluid flow. 

The PCM region was treated as an 


The effects of the non-steady inlet conditions Experimental 
of HTF on melting time, melting fraction, TES [105] 
capacity, solid-liquid interface, heat flux on 

tube surface and HTF outlet temperature 

were analyzed. 


The effect of natural convection of PCM 


during melting was taken into account 
with an effective thermal conductivity of 
the liquid phase of PCM. 


Adine and 2D 
Qarnia 
[106] 


Conservation energy 
equations 
Control volume approach 


developed. 


The flow of HTF was dynamically 


The axial conduction and viscous 
dissipation in the fluid were negligible. 
The thermophysical properties of the 


This parametric study could provide 
guidelines 

for system thermal performance and design 
optimization. 


Experimental 
[105] 


HTF and PCMs were independent of the 


temperature. 


The thermophysical properties of solid 


and liquid phases of PCM were equal, 
except their thermal conductivities. 


Apart from applying fins to PCMs, the effect of metal foams on 
heat transfer enhancement in rectangular encapsulated PCMs 
was investigated by Tian and Zhao [98] in 2011. The numerical 
investigation was based on the two-equation non-equilibrium 
heat transfer model, in which the coupled heat conduction and 
natural convection were considered at the phase transition and 
liquid zones. The numerical results were validated by experimen- 
tal data. The main findings of the investigation were that the heat 
conduction rate increases significantly by using metal foams, due 
to their high thermal conductivities, and that natural convection is 
suppressed owing to the large flow resistance in metal foams. 
In spite of this suppression caused by the metal foams, the overall 
heat transfer performance is improved when metal foams are 
embedded into PCM. This implies that the enhancement of heat 
conduction offsets or exceeds the natural convection loss. The 
results also indicated that for different metal foam samples, heat 
transfer rate can be further increased by using metal foams with 
smaller porosities and bigger pore densities. 


4.1.2. Cylindrical containment models 

Three modes of cylindrical PCM container configurations are 
classified. The first model is the PCM filling the shell and the heat 
transfer fluid (HTF) flowing through a single tube (pipe model). 
In the second mode the PCM fills the tube and the HTF flows parallel 
to the tube (tube-tube model). The third cylinder model is the PCM 
encapsulated in a shell and HTF flows through the tube outside the 
shell to improve the heat transfer of PCMs (tube-shell model). 

Ismail and Alves [119] numerically analyzed a pipe model with 
the PCM encapsulated in a shell while water as the HTF flowing 
through the tube. Radial conduction was assumed to be dominated 
during the process of solidification and the energy equation for the 


PCM was solved in conjunction with the equation governing 
the fluid bulk temperature as a function of time by employing 
the control-volume method. The effect of Biot number, ratio of 
diameters and inlet fluid temperature on the thermal performance 
of the storage unit were presented. Lacroix [105] developed an 
enthalpy based model to predict the transient thermal perfor- 
mance characteristics of a pipe LHTES with the PCM on the shell 
side and the HTF circulating inside the tube. Including the axial 
conduction effect in the PCM and employing an enthalpy based 
method, the coupled conservation equations governing the heat 
transfer in the PCM and HTF were solved iteratively using a finite- 
volume based finite-difference technique. Numerical results were 
presented for various thermal and geometric parameters such as 
shell radius, mass flow rate and inlet temperature of the HTF. 
It was concluded that for a given type of PCM, these parameters 
must be selected carefully in order to optimize the performance of 
the storage unit. Then the tube-tube model was theoretically 
studied by Esen and Durmus [120] in 1998. A theoretical model 
based on an enthalpy method was developed to predict the effects 
of various thermal and geometric parameters, configurations of 
the cylindrical container on the whole melting time for different 
PCMs. The theoretical results showed that the whole melting time 
of PCM depends on not only thermal and geometric parameters of 
the container, but also on the thermophysical properties of the 
PCM. It was found that a thicker PCM pack would result in longer 
melting time. 

During 2006, Regin et al. [121] experimentally and numerically 
studied the third model that the melting behavior of paraffin wax 
encapsulated in a cylindrical capsule. The numerical analysis 
carried out by using enthalpy method and the results were verified 
with the experimental data. The experiments have been done by 
the visualization technique without disturbing the actual process 
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of melting. Numerical results indicated that the melting process is 
chiefly governed by the magnitude of the Stefan number, phase 
change temperature range and the capsule radius. The analysis 
showed that the agreement between the analytical and experimental 
results is significantly improved when the results are obtained by 
considering the phase change temperature range and the natural 
convection in the liquid phase instead of only considering the 
conduction heat transfer. 


4.1.3. Spherical containment models 

Generally, the PCM in the sphere is subjected to constrained 
(the solid and liquid phases have the same density) and uncon- 
strained melting (the solid portion drops out of the melting layer 
due to its higher density). 

The constrained melting of the PCM within a spherical container 
was investigated by Khodadadi and Zhang in 2001 [4]. The 
computations were based on an iterative, finite-volume numerical 
procedure by using the primitive-dependent variables, whereby the 
time-dependent continuity, momentum and energy equations in 
the spherical coordinate system were solved. With the buoyancy- 
driven convection gradually dominated the heat transfer mechan- 
ism due to the growth of the melt zone, the solid PCM melted faster 
in comparison to the bottom zone. When the buoyancy effects 
became remarkable, as many as three time-dependent re-circulat- 
ing vortices were observed from the numerical simulation. The 
computational findings were verified through qualitative con- 
strained melting experiments using a high-Prandtl number wax. 

An experimental and computational investigation directed at 
understanding the role of buoyancy-driven convection during con- 
strained melting of the PCM inside a spherical capsule was reported by 
Tan et al. in 2009 again [122]. The computations were based on an 
iterative, finite-volume numerical procedure that incorporated a 
single-domain enthalpy formulation for simulation of the phase 
change phenomenon. A Darcy's Law-type porous media treatment 
linked the effect of the phase change on convection. The computa- 
tional predictions pointed to the strong thermal stratification in the 
upper half of the sphere that resulted from rising of the molten PCM 
along the inner surface of the sphere thus displaced the colder fluid. 
The measured temperature data and computational results near the 
bottom indicate the establishment of an unstable fluid layer that 
promotes chaotic fluctuations and is responsible for waviness of the 
bottom of the PCM. Regin et al. [123] also discussed an experimental 
and numerical study of a spherical capsule filled with a paraffin wax 
(melting temperature range of 52.9-61.6 °C) as the constrained PCM 
and placed in a convective environment during the melting process. 
Time-resolved temperature readings were obtained by using a rake of 
nine thermocouples that were placed on a horizontal plane passing 
through the center of the sphere. Keeping the initial temperature of 
the PCM constant, the temperature of the convective heat transfer 
fluid varied from 70 to 82°C corresponding to Stefan numbers of 
0.1143-0.2501. A 1D model for phase change by employing the 
enthalpy method and finite-difference formulation was developed. 
The model was then used to investigate the effect of capsule size and 
the Stefan number on the instantaneous heat flux, liquid fraction and 
melting time. 

Considering the spherical geometries, apart from several ana- 
lytical investigations of diffusion controlled phase change pro- 
cesses reported before 1980s, the first study on the unconstrained 
melting of the PCM within spherical container was carried out by 
Moore and Bayazitoglu [124] in 1982. A mathematical model was 
developed by assuming that the liquid remains at the fusion 
temperature during the whole period and solved by employing 
the perturbation method. The solid-liquid interface positions and 
the temperature profiles were determined for various Stefan and 
Fourier numbers. The experimental data agreed well with the 


simulation results. During 1987, Roy and Sengupta [125] per- 
formed a similar study and found an analytical solution for the 
melting rate at the lower surface of the solid PCM in contact with 
the heated surface. Later, Roy and Sengupta [126] performed 
another study on the unconstrained melting process inside a 
sphere. Two zones of the melting were identified: a thin melting 
layer at the bottom and a thicker layer at the top of the sphere. 
They found that a significant amount of melting took place at the 
upper surface due to the significance of natural convection, 
whereas in previous research [124,125], the effect of natural 
convection on the melting process was neglected. 


4.2. Microencapsulated PCM 


Micro-capsulation is a kind of tiny capsules (around 0.1- 
1000 um in diameter) [127]. Micro-encapsulation has a higher 
heat transfer rates as compared to that of macro-encapsulation 
[128,129] due to a substantially higher surface area to volume 
ratio. Meanwhile, microencapsulation can tolerate the change in 
volume during phase change process [130]. Microencapsulation 
hence is a prominent technology for use of PCM in building 
materials. With the advent of PCM implemented in gypsum board, 
plaster, concrete and/or other wall covering materials, thermal 
storage can be part of the building structure even for light weight 
buildings. 


4.2.1. Incorporation with building components 

PCMs can be incorporated with gypsum board, plaster, concrete 
and/or other wall covering materials for the temperature control 
application purpose. 

In general, storage density and phase change temperatures are 
very important parameters, since they determine the storage 
system size and human thermal comfort. Further, accurate estima- 
tion of the PCM enthalpy variation in the working temperature 
range is essential for correct mathematical modeling of the storage 
system | 131]. 

The first research on the application of PCMs in buildings has 
been carried out by Telkes in 1975 [132]. Na2S04-10H20 with a 
melting temperature of 32 °C was incorporated with walls, parti- 
tions, ceilings and floors to serve as temperature regulators. 

Determination criteria for the optimal phase change temperature 
and storage density were provided by Zhou [133] as following: 


Tm,opt = Ty + Q/htstor (26) 
T, =tg TgttnTn/tattn (27) 
Dopt = tnh(Tmopt a Tn)/pL (28) 


where Tmopt is optimal phase change temperature and Dopt is 
optimal thickness of the wall; Q represents the heat absorbed by 
unit area of the room surface; T, is the average room temperature; 
tstor is the heat storage time; the subscripts n and d represent night 
and daytime respectively. 

Kuznik et al. [134] developed a one-dimensional model to 
optimize the thickness of phase change material only considering 
pure conduction. The simulation was based on a lightweight wall 
and interior/exterior temperature evolutions within a period of 
24h. The results showed that an optimal PCM thickness exists. 
With the optimal PCM thickness of 1 cm, the thermal inertia of the 
building was doubled and the PCM can reduce the temperature 
fluctuations inside rooms. To provide guidelines useful in selecting 
an optimal PCM for building wallboard, a numerical model was 
developed to determine the optimal melting temperature of PCM 
by Neeper [135]. The results indicated that the optimal melting 
temperature depends on the average room temperature. When the 
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melting temperature of PCM closed to the average room tempera- 
ture, the maximum diurnal energy storage was achieved. Xiao 
et al. [136] established a simplified theoretical model to optimize 
an interior PCM for energy storage in a lightweight room. In this 
paper energy conservation equations were presented to calculate 
the optimal phase change temperature and the total amount of 
latent heat capacity. The analytical results agreed with Neeper's 
finding that the optimal melting temperature depends on the 
average indoor air temperature. Further, and it also depends on 
the radiation absorbed by the PCM panels. 

A 1-D enthalpy model for analyzing the thermal performance 
of a shape-stabilized PCM floor was developed by Xu et al. [137]. 
By using the model, the influence of various parameters (melting 
temperature of PCM, thickness of PCM layer, latent heat of fusion, 
and thermal conductivity of PCM) on the thermal performance 
was analyzed by using a fully implicit the finite-difference method. 
It was found that for a given position or weather condition, the 
suitable melting temperature of PCM is roughly equal to the 
average indoor air temperature of sunny winter days, the heat of 
fusion and the thermal conductivity of PCM should be larger than 
120 kJ/kg and 0.5 W/(m K), respectively, and the thickness of shape 
stabilized PCM plate used under the floor should not be larger 
than 20 mm. Zhou et al. [138] extended the investigation of the 
influence of various parameters from floor to walls and ceiling by 
using a verified enthalpy model. The model was validated by 
experimental work, and simulated results agreed well with the 
measured data. The numerical results showed that for the present 
conditions, the optimal melting temperature is about 20°C and 
the heat of fusion should not be less than 90 kJ/kg, thin PCM plates 
with large areas are advantageous and the effect of PCM plates 
located at the inner surface of interior wall was superior to that of 
exterior wall (the south wall). A similar result to Xu et al. [137] was 
also obtained for the condition of the PCM wall and ceiling that the 
thermal conductivity should be larger than 0.5 W/(mK). 

Athienitis et al. [139] experimentally and numerically studied 
the application of PCM in building envelope components for 
thermal storage. A 1-D non-linear enthalpy model was developed 
to simulate the transient heat transfer process in the walls using 
explicit finite difference scheme. The numerical result showed that 
the utilization of the PCM gypsum board may reduce the max- 
imum room temperature by about 4°C during the day time and 
can reduce the heating load at night significantly. 

The numerical simulation of passive heating with PCMs was 
carried out by Chen et al. [140] by using a 1-D nonlinear 
mathematical model. The effective heat capacity method was used 
to simulate the PCM heat transfer problem by using the software 
MATLAB. The result indicated that applying proper PCM to the 
inner surface of the north wall in the ordinary room can not only 
enhance the indoor thermal comfort dramatically, but also 
increase the utilization rate of the solar radiation. Consequently 
the heating energy consuming is decreased and the goal of saving 
energy has been achieved. The energy-saving rate of heating 
season can get to 17% or higher. 


4.2.2. Microencapsulated phase change slurry (MPCS) 
Microencapsulated phase change slurry (MPCS) as a new 
technique has been getting more and more attention since this 
technique serves as both the heat transfer fluids and energy 
storage media. Consequently it can increase the thermal storage 
capacity and improve the energy efficiency of the thermal system. 
Due to the MPCS' increasingly important role, some review papers 
on the properties, heat transfer characteristics and application of 
MPCS have been published in recent years [127,141-143]. As a 
kind of suspension, the heat transfer and fluid flow characteristics 
of MPCS are different from those of the marco-encapsulated PCMs. 


In accordance to the review of Delgado et al. [142], the heat 
transfer characteristics of the MPCS can be classified into forced 
convection and natural convection. 


4.2.2.1. Laminar forced convection heat transfer. The MPCS used to 
enhance convective heat transfer is generally laminar due to the high 
viscosity. Charunyakorn et al. [144] firstly conducted a numerical 
simulation of MPCS flow in circular tubes under different boundary 
conditions. The energy equation was formulated by taking into 
consideration both the heat absorption (and release) during phase 
transition and solved using an implicit finite difference method. Their 
parametric study showed that the bulk Stefan number and volumetric 
particle concentration are the two dominant parameters. Their model 
also showed that the heat transfer enhancement ratio of the fluid 
could be 2-4 times of that of water. Zhang and Faghri [145] modified 
the model of Charunyakorn [144] to include the effects of the 
microcapsule walls, initial sub-cooling and the melting temperature 
range. A temperature transforming model was used to solve the 
melting in the microcapsule by instead of a quasi-steady model. 
The numerical results showed a very good agreement with the 
experimental results of Goel et al. [146]. Their numerical results 
suggested that the differences between the experimental results of 
Goel et al. and numerical predictions of Charunyakom could be 
decreased greatly if these additional factors are considered. The 
authors however did not give any correlation or other similar 
criteria that could be used for future design. 

Complicated source terms are applied to the theoretical models 
above to account for the phase change effects. In order to simplify 
the simulation process, Alisetti and Roy [147] developed a simpler 
effective specific heat model to study the heat transfer in MPCS. The 
results showed that the exact form of the specific heat function is not 
critical as long as the latent heat is incorporated correctly within a 
finite melting temperature range. Based on the concept of effective 
specific heat, Roy and Avanic [148] developed a model to analyze the 
laminar forced convection heat transfer to a PCM suspension in a 
circular duct under constant wall heat flux condition. The model has 
been verified with experimental data. The results demonstrated that 
Stefan number is the only parameter that has a significant impact on 
the thermal performance. Sub-cooling effect only can affect the heat 
transfer performance at very low heat fluxes or when the inlet 
temperature is much lower than the melting point. A simple 
correlation to predict the wall temperature rise as a function of the 
tube length was also developed for future design. 

Hu and Zhang [149] introduced a model to provide a novel insight 
for the forced convective heat transfer enhancement of MPCS flowing 
through a circular tube with constant heat flux. The model was 
developed based on effective specific heat capacity method and used 
to analyze the influence of various factors on heat transfer enhance- 
ment. A modified Nusselt number was introduced since the conven- 
tional Nusselt number correlations for internal flow of single phase 
fluids are not suitable for accurately describing the heat transfer 
enhancement with MPCS. The modified local Nusselt number is 
inversely proportional to the dimensionless wall temperature. The 
numerical results indicated that the exact nature of the phase change 
process strongly affects the degree of convective heat transfer 
enhancement. The Stefan number and the PCM concentration 
resulted to be the most important parameters on the improvement 
of heat transfer in MPCS. Zeng et al. [150] analyzed experimentally 
and numerically the enhanced convective heat transfer mechanism 
of MPCS in the thermal fully developed range. An enthalpy model 
was proposed and developed for simulating the forced convective 
heat transfer characteristics of laminar flow with MPCS in a circular 
tube under constant wall heat flux. The results showed that in 
the phase change heat transfer region the Stefan number and the 
dimensionless phase change temperature range number are the 


672 S. Liu et al. / Renewable and Sustainable Energy Reviews 33 (2014) 659-674 


most important parameters influencing the Nusselt number fluctua- 
tion profile and the dimensionless wall temperature. The bulk fluid 
Reynolds number, particle diameter and the microcapsule volumetric 
concentration also influence the Nusselt number profile and the 
dimensionless wall temperature but are independent of the phase 
change process. 


4.2.2.2. Turbulent forced convection heat transfer. A majority of past 
studies have considered laminar forced convection heat transfer to 
MPCM and only a few numerical investigations on turbulent 
forced convection heat transfer have been reported. However, 
the turbulent flows occur in majority of engineering applications 
of MPCS. In order to investigate the turbulent heat transfer to PCM 
suspensions, Roy and Avanic [151] presented an effective specific 
heat capacity model for PCM suspensions in a circular tube with 
constant heat flux under turbulent flow condition. The numerical 
results were found to agree with previously published 
experimental data and showed that this effective specific heat 
model can effectively model turbulent heat transfer with PCM 
suspensions. The most primary parameter in heat transfer was the 
Stefan number. The melt temperature range and degree of sub- 
cooling are other two important parameters. Royon and Guiffant 
[152] developed a numerical model to describe the thermal 
behavior of a slurry of PCM slurries flow in a circular duct under 
different flow parameters (Reynolds number). The simulation 
results demonstrated that the influence of Reynolds number on 
the minimum length of the heat exchanger is relatively weak. 


4.2.2.3. Natural convection heat transfer. Inaba et al. [153] 
developed a 2-D model to study the fluid flow and heat-transfer 
characteristics for Rayleigh-Bénard natural convection of non- 
Newtonian PCM slurry. They found that Rayleigh number, 
Prandtl number and aspect ratio could be the primary factors for 
most of Newtonian and non-Newtonian fluids to evaluate the 
natural convection. A modified Stefan number was defined in the 
paper to evaluate the natural convection in a PCM slurry. Later, 
Inaba et al. [154] presented another 2-D model to study the 
natural convection heat transfer of a MPCS. It was found that the 
natural convection effect and heat transfer enhancement were due 
to the contribution of the latent heat transfer. 


4.3. Summary of numerical models 


It can be seen that the numerical models collected above for 
encapsulated PCMs all have their limitations. As the practical 
phase transformations in different applications are complicated, 
and the thermal conditions are not ideal. Various assumptions 
were set up for each numerical solution according to different 
numerical simulation purposes. There is none of numerical models 
can describe the phase change problem properly. Hence, the 
numerical modeling of LHTES only approximately approaches 
but not accurately presents the practical thermal performance of 
LHTES. Therefore, more investigation should be carried out to 
study the phase change problems so as to provide a more detailed 
insight of the actual heat transfer behaviors inside PCMs during 
the phase transformations. 


5. Conclusions 


This paper presented a comprehensive review of mathematical 
and numerical methods applied to the solutions of phase change 
problems. Firstly heat transfer mechanisms during the charging and 
discharging processes were discussed in this review. The heat 
transfer mechanisms play the most important role on the numerical 
results. Consequently, it is important to weight the percentages of 


conduction and convection heat transfer taken in each stage. Then, 
it presented the fundamental mathematical descriptions of the 
phase change phenomenon, the Stephan problem and Neumann 
problem. At last, a description of the numerical solutions for phase 
change problems based on considering only conduction and con- 
sidering also convection. The PCMs have been applied in thermal 
energy storage applications widely due to their attractive features. 
Various numerical models for encapsulated PCMs in terms of 
macro-encapsulated and micro-encapsulated PCM were also ana- 
lyzed in this paper. Nevertheless, the incorporation of PCMs in a 
particular application calls for an analysis that will enable the 
researcher to understand the heat transfer principle inside PCM 
correctly. From this survey, some further works are recommended: 


© Further investigation should be carried out to obtain adequate 
knowledge of the phase change problems so as to minimize 
the assumptions made in the numerical models. 

© Development of new numerical model to broaden its applica- 
tion in phase change problems and improve the accuracy of its 
results. 

© Assessment of the stability and convergence on the numerical 
results is always required, and the numerical predictions should 
always be validated by using appropriate experimental data. 

© Unified dimensionless parameters needs to be developed to 
make the gained knowledge applicable to other cases. 

© Life cycle analysis, both economical and energetic feasibility of 
PCM application should be performed as LHTES is a more 
expensive form of thermal storage, but few scientific publica- 
tions covered this matter. 
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